#import library
library(hydroGOF)
library(clusterSim)
library(sqldf)


inputPath = "C:/in2"
outputPath = "C:/out2/"
inputFile = "_merge.csv"

model_result = data.frame()
model_coeff = data.frame()

listFile = list.files(path = inputPath,pattern=inputFile,full.names = TRUE,recursive=TRUE)
Size = length(listFile)

for(i in 1:Size){
	
	filename = listFile[i]
	
	start_index = regexpr("_[1-9]",filename)+1
	end_index = regexpr("_merge.csv",filename)-1
	rt = substr(filename,start_index,end_index)
	
	data = read.csv(filename)
	data = subset(data,data$station!=3)
	
	cook_model = lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data)
	cookValue = cooks.distance(cook_model)
	data["cookValue"] = cookValue
	
	nSample = nrow(data)	
	nPredict = length(cook_model$coefficients)-1
	cutoff = 4/(nSample-nPredict-1)
	
	data_4np = subset(data,data$cookValue<=cutoff)
	model_4np = lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data_4np)
	data_4np["predict"] = model_4np$fitted.values
	
	attach(data_4np)
	np_samples = nrow(data_4np)
	np_re = sum(abs(predict-pm25)/pm25)/np_samples*100
	np_rmse = rmse(pm25,predict)
	np_r = cor(pm25,predict)
	
	np_mfb = (2*sum((predict-pm25)/(pm25+predict)*100))/np_samples
	np_mfe = (2*sum(abs(predict-pm25)/(predict+pm25)*100))/np_samples
	
	Samples = np_samples
	R2 = round(np_r*np_r,3)
	RMSE = round(np_rmse,3)
	RE = round(np_re,3)
	MFB = round(np_mfb,3)
	MFE = round(np_mfe,3)
	
	intercept_coff = model_4np$coefficients[1]
	aod_coff = model_4np$coefficients[2]
	temp_coff = model_4np$coefficients[3]
	avg_temp_coff = model_4np$coefficients[4]
	avg_rh_coff = model_4np$coefficients[5]
	avg_preci_coff = model_4np$coefficients[6]
	
	model_result = rbind(model_result,data.frame(rt,Samples,R2,RMSE,RE,MFB,MFE))
	model_coeff = rbind(model_coeff,data.frame(rt,intercept_coff,aod_coff,temp_coff,avg_temp_coff,avg_rh_coff,avg_preci_coff))
	
	data_4np_out_file = paste(outputPath,rt,"_data_4np.csv",sep="")
	
	write.csv(data_4np,data_4np_out_file)
	
	
}
write.csv(model_result,"C:/out2/npp_model_rs.csv")
write.csv(model_coeff,"C:/out2/npp_model_coeff.csv")



print("Done")